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Existing approximate Riemann solvers do not perform well when the grid is not aligned with strong 
shocks in the flow field. Three new approximate Riemann algorithms are investigated to improve solution 
accuracy and stability in the vicinity of strong shocks. The new algorithms are compared to the existing 
upwind algorithms in OVERFLOW 2.1. The new algorithms use a multidimensional pressure gradient based 
switch to transition to a more numerically dissipative algorithm in the vicinity of strong shocks. One new 
algorithm also attempts to artificially thicken captured shocks in order to alleviate the errors in the solution 
introduced by “stair-stepping” of the shock resulting from the approximate Riemann solver. This algorithm 
performed well for all the example cases and produced results that were almost insensitive to the alignment of 
the grid and the shock. 


Nomenclature 


a = speed of sound 

a = Roe averaged speed of sound 

e t = total energy per unit mass 

F = inviscid flux vector 

k p = pressure switch sensor 

k T = temperature switch sensor 

p = pressure 

S - wave speed 

t = time 

U = conserved variable vector 

u = velocity in the x direction 

u = contact discontinuity speed 

u = Roe averaged velocity in the x direction 

jc = spatial direction 

P = switching function 

X = Roe averaged eigenvalue of the inviscid flux jacobians 

p = density 

= computational coordinates 
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Subscripts 

spatial index 
left state 

contact discontinuity 
right state 


i 

L 

M 

R 


Superscripts 

n = temporal index 

* = intermediate state 


I. Introduction 

T he accurate computation of supersonic and hypersonic flows place conflicting demands on the formulation of 
inviscid numerical flux functions. A high speed numerical scheme must possess sufficient dissipation to capture 
strong shocks without developing overshoots and oscillations in the vicinity of the discontinuity. The scheme 
must also possess numerical dissipation that is much smaller than the physical viscosity to accurately compute 
boundary layers. The simultaneous satisfaction of these two requirements makes the computation of viscous 
hypersonic flows extremely challenging, particularly at the higher Mach numbers where many of the best schemes 
experience numerical difficulties. 

Most numerical algorithms developed for high speed flows are based on an approximate solution of the Riemann 
problem. Approximate Riemann solvers such as those developed by Roe 1 or Harten, Lax, and van Leer (HLL) 2 * are 
widely used for high speed flows because of their numerical accuracy and robustness. These approximate solvers 
are derived for one dimensional flow, and extensions to two and three dimensions introduce a possible source for 
numerical errors. Approximate Riemann schemes are highly sensitive to how well the shock is aligned with the 
grid. For cases in which the grid is not aligned with the shock, cross coupling between the Riemann problems in the 
different directions occurs. The strong jump in the conserved variables across the shock is “interpreted” by the 
Riemann solvers in the non-aligned directions as strong jumps in the other characteristic variables. This introduces 
numerical error for schemes that possess numerical dissipation that scales directly with the cell interface normal 
velocity in directions tangential to the shock where this component of velocity is low. The lack of dissipation in 
these tangential directions causes the cross-coupling perturbations to be either amplified (for schemes based on the 
use of anti-diffusion) or to be neutrally stable, i.e., the disturbances propagate without damping. In both cases the 
solution quality is severely degraded, resulting in non-physical features often referred to as “carbuncles”. 

Ad hoc fixes 2,3,4,5 to the eigenvalues based on satisfying entropy conditions have been proposed as a remedy to 
the difficulties with non-grid aligned shocks in the past, but all such approaches are problematic in satisfying the 
conflicting requirements for hypersonic shock capturing. Some of the fixes severely degrade the solution in the 
boundary layer. In addition, these entropy corrections merely attempt to add dissipation rather than directly 
manipulating the anti-diffusion terms. It is reasonable to expect that this effect would be especially aggravated in the 
case of tetrahedral grids 6 where the cell interfaces would likely be badly misaligned with the shock discontinuity. 
In this effort three new flux formulations are implemented into the OVERFLOW 2.1 flow solver and their ability to 
accurately compute flows with strong shocks and their efficacy is demonstrated on a variety of test cases. 

II. Theory 

The Godunov family of numerical schemes may be examined using the 1-D Euler equations 


dU dF „ 
dt dx 


( 1 ) 


Here U is the vector of conserved variables \p , pu , pe t ] T and F is the inviscid flux vector in the x direction given 
by 
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F = 


' P i 

pu 1 + p 
ipe t +p)u\ 


( 2 ) 


Eq. (1) can be discretized to obtain 


U? +1 

Ax 


(3) 


where At and Ax are temporal and spatial increments respectively. 

Roe Scheme 

The Roe 1 approximate Riemann solver defines the interface numerical flux as 


F Roe ~ 





where A is the Roe averaged inviscid flux Jacobian matrix given by 

a=LXr 


(4) 


(5) 


\ 2,3 = [w, u + a, u - a] are the Roe-averaged eigenvalues of the inviscid flux Jacobians and L and R are the left and 
right Roe-averaged eigenvectors of the inviscid fluxes. The Roe scheme has a well known deficiency in enforcing 
the entropy condition (i.e. in preventing nonphysical expansion shocks) when the eigenvalues of A vanish in regions 
near sonic points or stagnation points. This deficiency has traditionally been addressed by limiting the eigenvalues 
such that they do not fall below an arbitrary minimum value. The entropy fix implemented in OVERFLOW 7 uses 
the following function to select the eigenvalues used in the flux reconstruction. 



HLLEM Scheme 

Harten, Lax and Van Leer 2 proposed a family of Godunov schemes know as the HLL approximate Riemann 
solvers. The HLL formulation simplifies the solution to the Riemann problem at the cell interfaces by using 
averaged values for the interior states. The simplest HLL scheme assumes only one interior state between the fastest 
left and right running waves (S L and S R ). These waves represent the smallest and largest signal speeds with which 
information of the initial discontinuity is transported. A schematic (x-t plot) of the HLL formulation is depicted in 
Fig. 1. For Einfeldt’s version 8 of the HLL scheme (denoted as HLLE), the single state Riemann solution is given by 


UhllA 

i 


u L ifs L > 0 
U*ifS L <0<S R 
U R if S R < 0 


Here the intermediate state is given by 


S r U r -S l U l -(Fr-Fl) 


(V) 


( 8 ) 
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The flux for the HLLE scheme can be written as follows 




s + -s 


s + -s~ 


( 9 ) 


where 


S + = max^ ,0} 
S =min{5 i ,0} 


( 10 ) 


Einfeldt estimated the wave speeds for the fastest left and right running waves ( S L and S R ) as 

S R = maxjij ,(u + a) R } ^ 

S L =mm{i i ,(u- a ) L ) 

where \ =[w,w + a,w-a]are the Roe-averaged eigenvalues of the inviscid flux Jacobians. This scheme is robust 

but overly diffusive. In addition, all information about the contact discontinuity is lost. In order to rectify this, 
Einfeldt 9 developed the HLLEM scheme. The HLLEM scheme differs from the HLLE scheme by the addition of 
anti-diffusion terms to the linear characteristic fields as follows 


r HLLEM 


(u l ,u r )= 


s + f(u l )-s~f(u r ) | 

S + -S- 



-U L -^S p d p R p 


( 12 ) 


where tt p are the coefficients associated with the linear degenerated and nonlinear fields for characteristic waves 
which expand U R — U L in terms of the right eigenvectors of the Roe averaged flux Jacobians R p . 

d p =i p {0 R -U L ) (13) 

where L p is the left eigenvectors of the Roe averaged flux Jacobian. The anti-diffusion coefficients 8 P in Eq. (13) 
are defined to take out excess numerical dissipation in the linear degenerated fields 

s l =j=r—z 

I«1 + a (14) 

S 2 = S 3 = 0 


The value of U in Eq. (14) is an approximation for the speed of the contact discontinuity. The HLLEM scheme has 
numerical dissipation of the same order as the Roe scheme. This scheme can also be implemented without using 
logic statements. 

HLLC Scheme 

Toro 10 developed a two-state approximate Riemann solution that includes the contact discontinuity in the 
intermediate state. A schematic (x-t plot) of the HLLC two-state formulation is depicted in Fig. 2. The flux for the 
HLLC scheme can be written as 
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Frllc 


F L ifS L > 0 
} F*ifS L <0<S M 
'F* R ifS M <0<S R 
[ F R ifS R < 0 


The contact wave speed S M is given by 1 


S„ = 


_ Pr u r ($r ~ u r)~ Pl u l (&L ~ u l)+ Pl - Pr 


Pr{Sr ~ u r)~ Pl{$l ~ u l) 


(15) 


(16) 


The HLLC scheme will solve both isolated shock and contact waves exactly. The choice of wave speeds also 
guarantees that the method will be positively conservative 12 . Solutions with this scheme have been found to be 
comparable to an exact Riemann solver, only with a stronger enforcement of the entropy condition. This scheme 
does not allow direct access to the anti-diffusion terms. Implementation of the scheme does require logical 
statements. 


HLLE+ Scheme 

Park et al. 13,14 noted that the form of the HLLEM flux shown in Eq. (12) can be used to specify the HLLE, 
HLLEM, and Roe fluxes with the appropriate choice of wave speeds S R and S L in Eq. (11) and the contact 
discontinuity speed U in Eq. (14). The schemes can be defined using: 


HLLE: 


HLLEM: 


Roe: 


S R = max|i 2 ,(u + a) R j 
S L = minjijju - a) L j 
\u\ = QO 


S R = max[i 2 , (u + a) R j 


J _ 

Sr +$l 

' 

2 


Sr - A 2 

S L = ^3 


(17) 


(18) 


(19) 


Park and Kwan 13 proposed to combine the wave speed estimates from the HLLEM scheme with the contact 
velocity estimate from Roe’s scheme to form the HLLE+ flux function: 
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HLLE+: 


S R = max| A 2 ,( u + u )r j 

S L = min{i 3 ,(M - a) L } (20) 

u\ = | u 

Tramel and Nichols 12 and Park et al. 13 ’ 14 both noted that the anti-diffusion term 5 1 is responsible for many of the 
shortcomings which afflict the HLLE/Roe family of flux functions and used directional pressure sensors to reduce 
the magnitude of the anti-diffusion term in the presence of strong discontinuities. The directional pressure sensor is 
calculated along grid lines. The HLLEM scheme or Roe’s scheme can be made to transition to the HLLE scheme as 
the anti-diffusion is reduced. Since the switch is only activated in the vicinity of strong shocks, the switching does 
not affect other areas of the flow such as boundary layers. This allows hypersonic shocks to be captured without 
encountering the “carbuncle” phenomena occurring and without corrupting the low numerical dissipation of the 
scheme in boundary layers. The switch is implemented by defining the contact discontinuity speed U as 

|« | =/W|w|+0- P Park )« ( 2 1 ) 

where the switching function fi Park is defined by 


J\.otfk p ia p 

Pork | o q 0 tfo erw j se 


( 22 ) 


and the pressure switch sensor k p is defined as 


/ 

= ABS - 

v 


Pm ~ 2 Pi + Pm 
Pm + 2 Pi + Pm , 


(23) 


The pressure switch sensor is defined as the average of the sensors in each individual direction when applied to two- 
or three-dimensional flows. Ref. 13 recommends a value of 0.1 for the threshold value A p . 

All of the HLLE/Roe family of schemes can be cast in traditional artificial dissipation form given in Eq. (4) with 
the following choice of eigenvalues: 

*HLLE, i = + S~ )* u -2.0*(l-<5 1 )* (s + s- )y(s + - S ~ ) 

w ,2 = |s + + s~ )* (u + a)- 2.0 *(l-£ 2 )* (s + s- ))/(s + - 5- ) (24) 

A HLLE, 3 = ((s + + S - )* (u - a) - 2.0 * (l - £ 3 y (,S + 5- ))/(5 + - S - ) 


The linear eigenvalues in the Roe dissipation matrix transition from 


1*1 l /*2 * 2 \ 

\u\ —> — l u + a I 
11 2 a K ; 


(25) 


as j3 Park approaches zero. Alternately, when the Roe-averaged velocity u goes to zero the switching function 
reduces the anti-diffusion from one to 0.5 as j5 Park approaches zero. This modified contact discontinuity velocity 
provides enough dissipation to reduce or eliminate “carbuncles” in the vicinity of strong shocks. 

Higher Order Spatial Reconstruction 

This effort focuses on extending and upgrading the Roe and HLLC Riemann solvers currently implemented in 
OVERFLOW 2.1 15 . OVERFLOW 2.1 uses a Monotone Upstream-centered Schemes for Conservation Laws 
(MUSCL) approach to obtain higher order estimates of the primitive variables for the Riemann inviscid flux 
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routines. The code currently includes five options for flux limiters for these MUSCL extrapolations: Koren 16,17 , 
minmod 18 , van Albada 19 , WENO 20,21 , and WENOM 21,22 . OVERFLOW 2.1 uses a switch to locally reduce the 
spatial order of accuracy for the extrapolation for the minmod, van Albada, WENO, and WENOM limiters in the 
code to first order when strong gradients are encountered. The switch is currently evaluated in each of the 
computational directions. The switch is defined as 

S Wi =(A s -\.0)*MAX{k p ,k T ) 

(26) 

S Wi =MIN[S Wi , 1.0) 

(27) 

fiau=h -MAx(s Wti ,S Wt ,S WM f 

(28) 


where controls the sensitivity of the switch and A s is defined to be greater or equal to one. The parameter k T is a 
sensor defined using Eq. (23) based on static temperature rather than pressure. The switching function p oid ranges 
from zero to one. Large gradients cause the switching function to reduce to zero locally. Eq. (28) is used to 
“thicken” the sensor by using neighboring points. The switching function multiplies the higher order contribution to 
the local variable reconstruction and causes to reduce to first order in space near large discontinuities in pressure or 
temperature. 


III. OVERFLOW Implementation 

The code has been modified to use a multidimensional switch rather than the current directional switch by taking 
the maximum value of the Eq. (23) pressure sensor in each computational direction. The code currently 

uses the NAMELIST parameter DELTA to control the sensor sensitivity. The value of the threshold level of the 
switch (A s ) used in the code is set to 1.0 /DELTA. Values of DELTA between 1 and 10 (A^ varying from 1 to 0.1) 
have been found to be adequate to eliminate carbuncles and errors induced when strong shocks are present that are 
not aligned with the computational grid system. The multidimensional switch is used to locally reduce the spatial 
order of accuracy of the minmod, van Albada, WENO, and WENOM flux limiters for the Navier-Stokes equations. 
The new switch is defined as 


[ 0.0 if S Wi t >1.0/ DELTA 
0New | l — tanhjlO q>p^j otherwise 

where 

<p p = MAX\S W . t / DELTA, 0 j 


and S w is the maximum of kp for all of the neighboring points in all directions 


S w = MAX, . , J+1 [k n ) 

w i t j,k \ Pi >m>n / 


=j-hj+l 
=jt-l,£+l 


(29) 


(30) 


(31) 


The temperature is no longer included in the switching function to eliminate the possibility of large temperature 
gradients near the wall causing the switching function to activate. 

The existing Roe scheme in the code has been modified to become the HLLE+ scheme described in the previous 
section. The HLLE+ scheme uses the new switching function (Eq. (29)) rather than the switching function of Park 
(Eq. (22)). The existing HLLC scheme was modified to use the new multidimensional switching function to locally 
transition to the HLLE scheme using the following flux definition 

Fhllc+ = P HLLC+F HLLC + (l - Phllc+ ) f hlle (32) 
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where 


Phllc+ ~ MAX (0 New*® -4) (33) 

This scheme is similar in spirit to that advocated by Quirk 23 * The modified HLLC scheme has been designated 
HLLC+. The HLLE++ scheme was derived from the HLLE+ scheme by replacing the HLLE+ eigenvalues with a 
blending of the HLLE+ and a dissipative version of the Roe eigenvalues. The new eigenvalues are defined as 

^ HLLE++ = ft New^HLLE + (l ~~ P New 'ft Roe @4) 


The Roe eigenvalues ( A Roe ) are defined as 


4*k)= 


^Roe, 1 = 

MAX( 


^Roe, 2 = 

h(|m + 

a\,e) 

^Roe, 3 = 

hQm— 

a\,e) 

defined as follows 

W 

if 

w> 

Z 2 +S 2 

2s 

- if 

w<, 


(35) 


(36) 


and s was chosen to be 0.3 * QuJ + a). 

The switching function of Park (Eq. (22)) was found to cause numerical convergence to hang for both the 
HLLC+ and HLLE+ algorithms because of its abrupt transition between states. The old switching function (Eq. 

(28) ) was found to switch too quickly to the more dissipative HLLE algorithm even for small values of the pressure 
sensor kp. The new switching function (Eq. (29)) eliminated both of these problems. A comparison of the new (Eq. 

(29) ), old (Eq. (28)), and Park’s (Eq. (22)) switching functions is shown in Fig. 3. DELTA values between five and 
ten with the new switch should provide similar results to the Park switch. 

The relative advantages and di sadvantages of the schemes discussed in this paper are summarized as follows: 










Roe 

S Captures stationary grid-aligned shocks exactly 

* Expansion shocks (entropy fix required) 

S Reduced dissipation for boundary layers, etc. 

S Exact capture of stationary contact discontinuities 

* Susceptible to carbuncles 
HLLE 

S Captures stationary grid-aligned shocks exactly 
S No expansion shocks (no entropy fix required) 

* Excessive dissipation for boundary layers, etc. 

x Does not exactly capture of stationary contact discontinuities 
S Not susceptible to carbuncles 
HLLEM 

S Captures stationary grid-aligned shocks exactly 
S No expansion shocks (no entropy fix required) 

S Reduced dissipation for boundary layers, etc. 

S Exact capture of stationary contact discontinuities 
x Susceptible to carbuncles 
HLLC 

S Captures stationary grid-aligned shocks exactly 
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S No expansion shocks (no entropy fix required) 

S Reduced dissipation for boundary layers, etc. 

S Exact capture of stationary contact discontinuities 

* Susceptible to carbuncles 

• HLLE+ 

S Captures stationary grid-aligned shocks exactly 
S No expansion shocks (no entropy fix required) 

S Reduced dissipation for boundary layers, etc. 

S Exact capture of stationary contact discontinuities 

* Reduced Susceptibility to carbuncles 

• HLLC+ 

S Captures stationary grid-aligned shocks exactly 
S No expansion shocks (no entropy fix required) 

S Reduced dissipation for boundary layers, etc. 

S Exact capture of stationary contact discontinuities 

* Reduced Susceptibility to carbuncles 

• HLLE++ 

S Captures stationary grid-aligned shocks over « 2 grid points 
S No expansion shocks (no entropy fix required) 

S Reduced dissipation for boundary layers, etc. 

S Exact capture of stationary contact discontinuities 
S Much Reduced Susceptibility to carbuncles 

In summary, all of the approximate Riemann algorithms discussed here have been modified to reduce their 
susceptibility to carbuncles. All of the algorithms use a switch to locally reduce the order of extrapolation to first 
order in the presence of strong gradients. The HLLC and Roe scheme use a directional switch based on pressure and 
temperature. The baseline HLLC has no further corrections. The original Roe algorithm also limits the eigenvalues 
in the anti-diffusion terms to impose an entropy fix to eliminate expansion shocks. The HLLE+ and HLLC+ 
algorithms use a multidimensional pressure gradient-based switching function to change to the more dissipative 
HLLE scheme near shocks. The HLLE++ is the FILLET scheme with the addition of a limiter invoked by the 
pressure gradient switch applied to the eigenvalues. 

TV. Results 

Solutions using the current OVERFLOW 2.1 HLLC and Roe schemes are denoted as HLLC and Roe in this 
section. Solutions using the new schemes are denoted as HLLC+, FILLET, and HLLEtt. Note that HLLC+ is 
exactly the same as HLLC when DELTA=\. Also the HLLEtt scheme is identical to the FILLET scheme when 
DELTA=\. All solutions were run with 3 rd order spatial operators unless otherwise noted. The 5 th order spatial 
WENO schemes behave in a similar manner to the 3 rd order algorithms. 

Inviscid Supersonic Ramp 

The first case demonstrates the problem with the current Riemann solvers in OVERFLOW 2.1 and shows the 
improvement that result from using the new schemes. The inviscid M=4.0 flow over a two-dimensional 30° ramp is 
examined. All solutions were performed using the van Albada flux limiter. Fig. 4 shows the 59x219 grid colored by 
density for the HLLETT scheme with DELTA=5. Density contours for the HLLC+ algorithm with DELTA=\ are 
shown in Fig. 5a. The shock is not aligned with the grid. The density contours at the shock are seen to “stair-step” 
across the grid and errors are produced that are convected downstream of the shock. Similar results were obtained 
with the HLLE+ and HLLE++ algorithms with DELTA = 1 and for the Roe and HLLC algorithms with all DELTA 
values. Density contours for the HLLC+ algorithm with DELTA=5 are shown in Fig. 5b. The error streaks behind 
the shock are significantly reduced. Similar results were obtained with the HLLE+ algorithm. Density contours for 
the HLLETT algorithm with DELTA=5 are shown in Fig. 5c. The error streaks are almost completely removed by 
this algorithm. 

The density on the comparison line indicated in Fig. 5 are shown for all five algorithms with DELTA=l 9 5, 10, 
and 20 in Fig. 6. All of the algorithms display at least some of the convected error that is originating from the shock 
when DELTA=\. Increasing DELTA has no effect on reducing the error with the Roe or HLLC algorithms. 

9 

American Institute of Aeronautics and Astronautics 



19 th AJAA Computational Fluid Dynamics Conference 
22-25 June 2009, San Antonio, TX 


AIAA-2009-3988 


Increasing DELTA is shown to greatly reduce the erroneous wiggles for the HLLE+ and HLLC+ algorithms. The 
HLLE-H- algorithm decreases the error even further when DELTA is increased. These results show that the post- 
shock oscillations are a result of the multi-dimensional resolution of the shock since the error is significantly 
reduced with the use of a multidirectional switch. 

It has often been conjectured that reducing the spatial order of accuracy could remove the error associated with 
the misaligned shock through the addition of large amounts of numerical dissipation. The density on the comparison 
line for all algorithms using first order spatial algorithms is shown in Fig. 7. The density oscillations are reduced 
compared to the 3 rd order spatial results for all algorithms, but the error is still much larger than that seen for the 3 rd 
order spatial HLLC+, HLLE+, and HLLE++ algorithms with DELTA>\. Hence simply reducing to a first order 
spatial algorithm at the shock will reduce but will not eliminate the error generated at the shock. 

The minimum oscillation defined in terms of the difference of the maximum and minimum density divided by 
the average density for each algorithm is shown in Table 1. The HLLE++ algorithm produces the lowest density 
oscillation. The average value of the density for the HLLC+, HLLE+, and HLLE++ algorithms are in reasonable 
agreement with the theoretical density behind the shock (p/p o0 =3.703). The HLLC+, HLLE+, and HLLE++ 
algorithms all produce results with lower density oscillations and better agreement with the theoretical density than 
the central difference algorithm. 


Algorithm 

DELTA 

pip*, 

AP^ave (%) 

Roe 

10 

3.682 

17.19 

HLLC 

10 

3.679 

22.77 

HLLC+ 

20 

3.684 

0.573 

HLLE+ 

20 

3.683 

0.706 

HLLE-H- 

20 

3.683 

0.360 


Table 1 . Minimum density oscillation for each algorithm for the M=4 30° ramp. 

The solutions were run using grid sequencing with 100 iterations on the coarse mesh, 100 iterations on the 
medium grid, and finally 1000 iterations on the fine grid. The convergence history for the Roe and HLLE++ 
algorithms is shown in Fig. 8. The HLLC convergence history is similar to the Roe algorithm. The HLLC+ and 
HLLE+ convergence history is similar to the HLLE++ results. The HLLE++ algorithm with DELTA>\ has much 
better convergence properties than the algorithm with DELTA =1. The Roe and HLLC algorithm convergence does 
not improve with increasing values of DELTA. 

The switching function J3 New is shown in Fig. 9 for the HLLE++ algorithm with DELTA= 5. The switching 
function rapidly approaches zero in the vicinity of the shock. The HLLC+ algorithm is equivalent to the HLLC 
algorithm in smooth regions of the flow, but the HLLC+ algorithm switches to the HLLE scheme near shocks where 
the switching function is zero. Thus it follows that the HLLE scheme is the agent for removing the error associated 
with the non-grid aligned shock for this case. 

Inviscid Supersonic Ramp with Channel 

The second case builds upon the first case with the addition of shock-shock interaction. The inviscid M=3.0 flow 
over a two-dimensional 30° ramp in a channel is examined. The channel height is twice the ramp height. All 
solutions were performed using the van Albada flux limiter. Fig. 10 shows the 95x101 grid colored by density for 
the HLLE-H- scheme with DELTA=5. Density contours for the HLLC+ algorithm with DELTA=\ are shown in Fig. 
11a. The shock is not aligned with the grid. The density contours are seen to possess streaks downstream of the 
oblique and normal shocks. Similar results were obtained with the HLLE+ and HLLE-H- algorithms with 
DEL TA= 1.0 and for the HLLC algorithm with all DELTA values. No steady state solution could be obtained for the 
Roe algorithm with DELTA=\. Density contours for the HLLC+ algorithm with DELTA=5 are shown in Fig. 1 lb. 
The density contours behind the oblique shock are greatly improved, but there are still nonphysical density 
oscillations behind the normal shock. Similar results were obtained with the HLLE+ algorithm. Density contours 
for the HLLE++ algorithm with DELTA=5 are shown in Fig. 11c. There is no visible density streaking or 
oscillations present behind either the normal or oblique shock. 

The density on the comparison line indicated in Fig. 11 are shown for all five algorithms with DELTA=\> 5, 10, 
and 20 in Fig. 12. All of the algorithms display at least some of the convected error that is originating from the 
shock when DELTA=\. Increasing DELTA has no effect on reducing the error with the Roe or HLLC algorithms. 
Increasing DELTA is shown to greatly reduce the erroneous wiggles for the HLLE+ and HLLC+ algorithms. The 
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HLLE++ algorithm eliminates the error when DELTA is increased. All of the algorithms show differences in the 
expansion region when DELTA is increased. 

The solutions were run using grid sequencing with 100 iterations on the coarse mesh, 100 iterations on the 
medium grid, and finally 1000 iterations on the fine grid. The convergence history for the HLLC and HLLE++ 
algorithms is shown in Fig. 13. The Roe convergence history is similar to the HLLC algorithm. The HLLC+ and 
HLLE+ convergence history is similar to the HLLE++ results. The HLLE++ algorithm with DELTA>\ has much 
better convergence properties than the algorithm with DELTA =1. The Roe and HLLC algorithm convergence 
improves some with increasing values of DELTA. 

The switching function J5 New is shown in Fig. 14 for the HLLE++ algorithm with DELTA=5. The switching 
function rapidly approaches zero in the vicinity of the shock. The switch approaches zero in the expansion region 
by the ramp upper comer. 

Shock Tube 

Shock tube problems have long been used as test cases for unsteady flows and are used here to demonstrate that 
the algorithms are time accurate when DELTA > 1 . The position of the shock, expansion, and contact discontinuity 
in inviscid flow can be predicted theoretically. The conditions chosen for this problem are: 


Left state: p/p Q0 = 1.0, p/p 00 = 1.0 
Right state: p/p x = 0.1, p/p^ 0.1 

The results are evaluated at the nondimensional time of 0.2 (320 iterations for the nondimensional time step of 
0.000625). The solution was run using 401 grid points with ten Newton subiterations per time step. The density at 
t=0.2 for each algorithm is shown in Fig. 15. The solutions are almost indistinguishable and are in good agreement 
with theory for all DELTA values. 


Inviscid Isentropic Vortex Convection 

The ability to conserve both the vortex shape and strength is important in many unsteady cases in which a shed 
vortex interacts with bodies well downstream of the vortex origin. The following test case 25 can be used to examine 
the relative level of numerical dissipation and dispersion for inviscid flux algorithms. An isentropic vortex of non- 
dimensional strength T= 5 is centered on a uniform grid of size 10 units by 10 units. A 101x101 grid (Ax=0.1) was 
used for most of the results presented here. The vortex is allowed to convect downstream at free stream Mach 


number of 0.5 with a non-dimensional time step 


f 

dt 

V 


u 




zf J 


of 0.01. 


The velocity, temperature, density, and 


pressure for the vortex are given by 


= - ^~{ z - z 0 )exp[o.5(l - R 2 )] 


2 n 


w - (x - x 0 )exp[o.5(l - R 2 )] 
2n 


T = T- 


%yn 


expl 


M J ) 


p = T y ~ x 
P = P r 


(34) 


where R 2 = (x - x 0 ) 2 + (z - z 0 ) 2 and x 0 and z 0 define the location of the vortex center at any given point in time. 

Here y is the ratio of specific heats. The grid is given periodic boundary conditions in the flow direction. The 
vortex should complete one cycle on the grid every 1000 time steps. All calculations used 2 nd order time and three 
Newton subiterations. The vortex is allowed to convect through the grid five times. The grid and the starting vortex 
can be seen in Fig. 16. 
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The effect of spatial order of accuracy on the minimum density of the vortex for the HLLC scheme is shown for 
reference in Fig. 17. Note the dramatic difference in dissipation level between 1 st and 2 nd order space for this case. 
The minimum density for the vortex with all of the flux schemes and varying DELTA is shown in Fig. 18. The Roe 
and HLLC schemes remain 3 rd order in space for all values of DELTA. The HLLC+, HLLE+, and HLLE++ 
algorithms are relatively insensitive to the value of DELTA , with both algorithms reducing to about second order for 
DELTA=20. The HLLC+, HLLE+, and HLLE++ algorithms remain near 3 rd order in the preferred range of 1< 
DELTA <10 shown to adequately remove the error associated with non-grid aligned shocks for the previous cases. 

Hypersonic Cylinder Bow Shock 

The hypersonic bow shock experiment data of Holden, et. al 26 provides a good test of the ability of the numerical 
schemes to capture extremely strong shocks and to predict heat transfer at high speeds. The experiment was run at a 
free stream Mach number of 16.01, a Reynolds number (based on cylinder diameter) of 9.11xl0 4 , a free stream 
temperature of 77.8° R, and a wall temperature of 540° R. The nitrogen test medium can still be considered a perfect 
gas at these conditions. The flow over the 3 in. diameter cylinder is laminar for this case. This case has a 
temperature ratio across the shock of about 50 and a pressure ratio across the shock of about 300. Two 
computational grids were used in this study. Both grids had 101x101 with a wall spacing corresponding to y + value 
of 0.1. No effort was made to align the grid with the shock for one of the grids. The bow shock in the nose region 
was aligned with the grid for the second mesh. 

The flow field was initialized with a free stream Mach number of 16.01. Two levels of grid sequencing were 
used to set up the flow field (50 iterations on the coarse level and 50 iterations on the mid level). The case was run 
time accurately with a time step of 2.9x1 0' 5 sec. Three Newton subiterations were used at each time step with all of 
the inviscid flux algorithms. The grid and temperature contours for both grids are shown in Fig. 19. 

Temperature contours using the HLLC+ algorithm with DELTA= 1 and 5, and the HLLE++ algorithm with 
DELTA= 5 are shown in Fig. 20. The contours on the centerline behind the shock are slightly different between the 
DELTA=\ and DELTA=5 results. There are further differences as the flow inside the shock proceeds around the 
cylinder. The HLLE++ algorithm with DELTA= 5 produces much smoother and rounder contours than the HLLC+ 
algorithm. The Roe, HLLC, and HLLE+ results are similar to the HLLC+ results. 

Predicted pressure coefficient is compared to modified Newtonian theory for the HLLE++ scheme in Fig. 21. 
All of the algorithms produced similar results and are in good agreement with theory in the nose region where the 
theory is valid. Predicted heat transfer is compared to data in Fig. 22 for the HLLE+ scheme. All of the flux 
algorithms produced similar results and were in good agreement with the data. The stagnation heat transfer is 
slightly higher when DELTA= 1 on the shock aligned grid for all the algorithms. The stagnation heat transfer is 
slightly lower when DELTA= 1 on the non-shock aligned grid. All of the results are well within the scatter of the 
data. The solutions for DELTA> 1 all predict about the same level of heat transfer. Since all of the algorithms show 
the similar behavior in peak heat transfer, this can be attributed to the reduction of order of the flux limiters in the 
vicinity of the shock and not to the switching to the HLLE scheme. 

The L2 norm of the residuals is shown in Fig. 23 for the HLLE-H- algorithm on both grids. Convergence is 
improved with increasing DELTA for this algorithm. The solution converges an order of magnitude further on the 
shock aligned grid. Similar behavior was noted for all algorithms. The switching function j3 New is shown in Fig. 24 
for DELTA=5. The switching function is one except near the shock and at the downstream boundaries. 

Hypersonic Double Cone 

The hypersonic double cone experiments of Holden and Wadhams 27 provides a very demanding test of the ability 
of the code to predict the heat transfer at high speeds. The geometry for this case is shown in Fig. 25. Details of the 
flow field are shown in Fig. 26. The attached shock from the first cone interacts strongly with the detached shock 
associated with the second cone. This shock/shock interaction produces a transmitted shock which in turn impinges 
on the surface of the second cone. This impingement produces very high surface heat transfer rates and pressures. In 
addition, the high pressures which result at the cone-cone junction cause the flow to separate in this region. This 
separation bubble in turn interacts with the inviscid flow field, which impacts the strength of the transmitted shock, 
etc. and the entire flow field is very dependent on the strengths of the relative interactions. Computational 
predictions of this type of flow are very sensitive to numerics, grid resolution, etc. (Druget, Candler, and 
Nompelis 28 ). 

The computational results are compared to the experimental results of Holden’s Run 35. The experiment was run 
at a free stream Mach number of 12.06, a Reynolds number of 6.79xl0 4 per foot, a free stream temperature of 182° R 
and a wall temperature of 533° R. The test medium was pure nitrogen and is considered a perfect gas at these 
conditions. The flow conditions were chosen such that the flow over the entire length of the cone remains laminar. 
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Computations were run on a 305x129 grid. This grid is one level coarser than the medium grid used in Ref. 1 5. The 
grid colored by temperature is shown in Fig. 27. The bow shock on the second cone is not aligned with the grid. 

Temperature contours are shown in Fig. 28. The HLLC+ DELTA= 1 and HLLC+ DELTA=5 results are similar. 
Both of these cases show some “stair-stepping” behind the bow shock on the second cone with error that is 
convected downstream of the shock. The HLLE++ contours are much cleaner in this region of the flow. The Roe, 
HLLC, and HLLE+ results are similar to the HLLC+ contours. 

The surface pressure coefficient is shown for all the algorithms in Fig. 29. The surface heat transfer is shown in 
Fig. 30. The DELTA=1 results compare the best with the data for all algorithms. The Roe and HLLC results are 
very sensitive to the value of DELTA. The HLLC+, HLLE+, and HLLE++ results are less sensitive to DELTA for 
DELTA <10 than the Roe and HLLC algorithms, but a dependency still exists. The separation region at the 
intersection of the two cones decreases as DELTA is increased causing the surface pressure coefficient and heat 
transfer predictions to move away from the data. The final expansion region following the end of the second cone 
shows a large dependency on DELTA for all algorithms. This dependency can be attributed to the switching 
function and indicates that using values of DELTA between one and five would be preferable for this case. 

The switching function J5 New is shown in Fig. 31 for DELTA=5. The switching function is less than one in the 
vicinity of the shock, at the upper comer at the start of the expansion region, and at the start and end of the 
separation region at the intersection of the two cones. Convergence for the Roe and HLLE++ algorithms are shown 
in Fig. 32. Both of these algorithms indicate convergence plateaus with a higher error level for DELTA=20. 

Capsule 

The final case is a three dimensional capsule geometry shown in Fig. 33. The flight conditions were a Mach 
number of 4.02, 16° angle-of-attack, a free stream temperature of 342° R, and a Reynolds number of 2.49x1 0 5 per 
inch. The ratio of specific heats was 1.246. The walls were treated using adiabatic no-slip boundary conditions. 
Two grid systems with 7.3x1 0 6 points were used in this study and are shown in Fig. 34. One grid was aligned with 
the bow shock while the other grid was not aligned with the shock. The non-shock aligned grid is representative of 
the grid required for a moving body simulation where the shock will move during the simulation. 

Mach number contours on the symmetry plane for the shock aligned grid are shown in Fig. 35 for the HLLC+ 
scheme with DELTA=\ , HLLC+ scheme with DELTA= 5 and the HLLE++ scheme with DELTA=5. The HLLC+ 
scheme with DELTA=\ shows carbuncle structures along the bow shock. All of the schemes had similar carbuncles 
when DELTA=\. The Roe and HLLC algorithms produced carbuncles for all values of DELTA. The HLLC+, 
HLLE+, and HLLE++ algorithms all produced reasonable results for DELTA>\. Mach number contours on the 
symmetry plane for the non-shock aligned grid are shown in Fig. 36 for the HLLC+ scheme with DELTA=l 9 
HLLC+ scheme with DELTA=5 and the HLLE++ scheme with DELTA=5. The HLLC+ scheme with DELTA=l 
shows carbuncle structures along the bow shock. All of the schemes had similar carbuncles when DELTA=\. The 
Roe and HLLC algorithms produced carbuncles for all values of DELTA. The HLLC+ and HLLE+ algorithms have 
a small amount of error from the shock “stair-stepping” in the region behind the bow shock for DELTA>\. The 
HLLE++ algorithm predicts a flow field that is almost identical to the shock-aligned grid results for DELTA>\. 

Pressure coefficient distributions on the symmetry plane are shown in Fig. 37 for the shock-aligned grid and in 
Fig. 38 for the non-shock aligned grid. The distributions show the same trends as seen in the Mach number 
contours. A rise in pressure is seen at about Z=-0.6 on the windward side for the non-shock aligned grid. This 
pressure rise is associated with a disturbance in the bow shock. Fig. 39 shows the pressure coefficient on the 
symmetry plane for the HLLC+ algorithm with DELTA=5 and the HLLE++ algorithm with DELTA=5 on both grids. 
All of the solutions are in agreement except for the HLLC+ result with DELTA=5 on the non-shock aligned grid. 

The switching function jS New is shown in Fig. 40 for DELTA=5 for both grids. The switching function behaves 
similarly for both grids. The switching function is less than one in the bow shock region and in the vicinity of the 
weaker shocks on the leeward side of the capsule 


V. Conclusions 

Traditional approximate Riemann solvers such as the HLLC or Roe algorithms in OVERFLOW 2.1 do not 
perform well when the grid is not aligned with strong shocks in the flow field. Numerical error in the form of 
“carbuncles” or “stair-steps” is formed in the vicinity of the shock. This error is convected downstream and reduces 
the accuracy and stability of the solution. In high speed moving body simulations and complex flow fields it is often 
impossible to modify the grid to maintain shock alignment. 

Four approaches for removing the non-grid aligned shock error. These approaches include adding entropy fixes, 
locally reducing to first order in space in the shock region, modifying the anti-diffusion terms of the flux, and locally 
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switching to a more dissipative algorithm near the shock. None of these fixes by themselves were capable of 
eliminating the error introduced by non-grid aligned shocks. The HLLE++ scheme includes all of the approaches 
except for the entropy fix and was found to perform better than any of the other algorithms investigated. 

The HLLC+, LILLEE, and HLLE++ algorithms were shown to provide significantly better solutions than the 
original Roe or HLLC algorithms for problems with strong shocks that are not aligned with the grid. This is due to 
the fact that the HLLC+, HLLEE, and HLLE++ use a pressure gradient based switching function to locally transition 
to the more dissipative HLLE algorithm in the vicinity of strong shocks. The results here also indicate that the use 
of a temperature sensor in addition to the pressure sensor is not required to obtain stable and accurate solutions with 
large flow field discontinuities. The temperature sensor may be invoked in the boundary layer, while the pressure 
sensor would not be invoked near the wall. 

The HLLC+, HLLEE, and HLLE++ were also shown to maintain time accuracy and provide accurate heat 
transfer. The new algorithms were also shown to provide low numerical dissipation for the convecting vortex case 
for the range of DELTA values required to remove error in the vicinity of non-grid aligned shocks. The HLLE++ 
algorithm produced almost grid independent solutions regardless of grid alignment with strong shocks. Numerical 
convergence was improved for all of the new algorithms for \>DELTA<\0. 

The switching function was shown to affect expansion regions as well as shocks in the example cases. This is 
not a desirable feature. Fortunately low values of the DELTA input parameter (\<DELTA<5) were found to be 
sufficient to minimize the formation of the numerical error near the shock for this algorithm. The lower values of 
DELTA will also assure that the algorithm maintains low numerical dissipation. Hence it is recommended that 
smaller values of DELTA be used in applications with the HLLC+, HLLEE, or HLLE-H- algorithms. 

The additional eigenvalue limiting used in the HLLE++ scheme consistently produced the best results for all the 
problems investigated here. The eigenvalue limiting in the HLLEef scheme tends to smear the shock over two 
cells. The HLLC+ and HLLE+ schemes are capable of capturing the shock in one cell. Based on the results here, 
the accuracy and stability of the results more than offset the slight smearing of the shock with the HLLE++ scheme. 

All of the examples in this paper used the 3 rd order spatial MUSCL extrapolation to reach higher order. The 5 th 
order WENO extrapolations used in OVERFLOW 2.1 have similar sensitivity to carbuncles and “stair-stepping” as 
the 3 rd order spatial MUSCL schemes. The WENO schemes show the same improvements when used in 
conjunction with the HLLC+, HLLEE, or HLLEee algorithms. 
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Figure 3. Switching functions for the HLLC+ and HLLE+ algorithms. 
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Figure 4. Grid colored by density for the M=4.0 30° ramp. 



Figure 5a. Density contours using HLLC+ with DELTA=\ for the M=4.0 30° ramp. 
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Figure 5b. Density contours using HLLC+ with DELTA = 5 for the M=4.0 30° ramp. 



Figure 5c. Density contours using HLLE++ with DELTA=5 for the M=4.0 30° ramp. 
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Figure 6. Density distribution on the comparison line for the M=4.0 30° ramp. 
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Figure 7. Density distribution on the comparison line using 1 st order spatial fluxes for the M=4.0 30° ramp. 




Figure 8. Convergence history for the M=4.0 30° ramp. 
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Figure 9. j3 New for the M=4.0 30° ramp for DELTAS. 



Figure 10. Grid colored by density for the M=3.0 30° ramp in a channel. 
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Figure 11a. Density contours using HLLC+ with DELTA=\ for the M=3.0 30° ramp in a channel. 



Figure 1 lb. Density contours using HLLC+ with DELTA = 5 for the M=3.0 30° ramp in a channel. 



Figure 11c. Density contours using HLLE++ with DELTA = 5 for the M=3.0 30° ramp in a channel. 
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Figure 12. Density distribution on the comparison line for the M=3.0 30° ramp in a channel. 
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Figure 13. Convergence history for the M=3.0 30° ramp in a channel. 



Figure 14. fi New for the M=3.0 30° ramp in a channel for DELTA=5. 
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Figure 15. Density distribution in a shock tube. 
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Figure 16. Initial density contours for the convecting vortex. 
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Figure 17. Effect of spatial order for the inviscid flux terms for the converting vortex. 
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Figure 18. Minimum density for the convecting vortex. 



Figure 19. Grids colored by temperature. 
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Figure 20. Temperature distribution for the M=16.01 cylinder bow shock. 
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Figure 21. Pressure coefficient for the M=16.01 cylinder bow shock. 
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Figure 22. Heat transfer for the M=16.01 cylinder bow shock. 




Figure 23. Convergence history for the M= 16.01 cylinder bow shock. 
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Figure 24. j3 New for the M= 16.01 cylinder bow shock for DEL TA =5. 



Figure 25. Geometry for the hypersonic double cone (all dimensions in inches). 
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Figure 27. Grid colored by temperature for the hypersonic double cone. 
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Figure 28. Temperature contours for the hypersonic double cone. 
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Figure 29. Pressure coefficient for the hypersonic double cone. 
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Figure 30. Heat transfer for the hypersonic double cone. 



Figure 31. fi New for the hypersonic double cone for DELTAS. 
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Figure 32. Convergence history for the hypersonic double cone. 



Figure 34. Grid on symmetry plane for capsule. 
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Figure 35. Mach number on the symmetry plane for the shock aligned grid. 



Figure 36. Mach number on the symmetry plane for the non-shock aligned grid. 
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Figure 37. Pressure coefficient on the symmetry plane for the shock aligned grid. 
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Figure 38. Pressure coefficient on the symmetry plane for the non-shock aligned grid. 



Figure 39. Pressure coefficient on the symmetry plane for both the shock aligned and non-shock aligned grids. 
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Figure 40. 0 New for the capsule for DELTA=5. 
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